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Abstract 

A constructive approach to theory of diffusion processes is proposed, which 
is based on apphcation of both the symmetry analysis and method of mod- 
elhng functions. An algorithm for construction of the modelling functions is 
suggested. This algorithm is based on the error functions expansion (ERFEX) 
of experimental concentration profiles. The high-accuracy analytical descrip- 
tion of the profiles provided by ERFEX approximation allows a convenient 
extraction of the concentration dependence of diffusivity from experimental 
data and prediction of the diffusion process. Our analysis is exemplified by its 
employment to experimental results obtained for surface diffusion of lithium 
on the molybdenum (112) surface pre-covered with dysprosium. The ERFEX 
approximation can be directly extended to many other diffusion systems. 



1 Introduction 



Experimental and theoretical studies of diffusion processes are of a great importance 
for various branches of physics, biology, chemistry and other natural sciences. In ad- 
dition, such studies have important applications in medicine and many technological 
processes. A special interest is exited by surface diffusion processes which appear in 
many physical and chemical systems. In particular, they are used in various kinds of 
nanotechnologies. 

The theory of diffusion processes started in 1855 when Pick derived his classical 
diffusion equation [Ij 



which still is a corner stone of the diffusion theory. In equation (fT]) is a diffusion 
coefficient, in general case depending on species concentration 9, and Xa with a = 
1, 2, 3 are spatial variables (summation over the repeated indices a is imposed). Being 
supplemented by an appropriate initial data, equation ([T]) serves as a background for 
description of such diffusion processes which are characterized by diffusion flows linear 
in concentration gradients and not depending explicitly on space and time variables. 
Two standard problems of a diffusion theory are: 

1) To describe time evolution of the diffusion process, and 

2) To specify the dependence of the diffusion coefficient on concentrations of dif- 
fusing species. 

Of course, these problems are closely related, since if we know how the diffusion 
coefficient depends on concentration 6, then the time evolution of the corresponding 
diffusion process can be found using the Fick equation ([T]) and the related initial data. 
On the other hand, if we know 6' as a function of time variable t and spatial variables 
Xa, then we can find D solving the inverse diffusion problem using again equation ([T]). 

Both mentioned problems are very complicated and in general need rather sophis- 
ticated techniques. Even if we know the diffusion coefficient as an explicit function of 
concentration, then generally speaking it is possible to find only an approximate (nu- 
merical) solution of the first problem if at all. The second problem has a much more 
complex character, but in the case of a sharp step-like initial 6 profile it is possible 
to use the Boltzmann-Matano (BM) approach [2] and reconstruct the concentration 
dependence D{6) of the diffusion coefficient. This approach enables one to make a 
numerical calculation of the diffusion coefficient, but its accuracy is not very high, 
especially for small and large concentrations 6. 

Experimental data and numerical solutions are very important for description of a 
diffusion process, but to formulate its theory it is desirable to create some analytical 
expressions for studied values. Unfortunately, there are only few known exactly solv- 
able realistic diffusion problems, the most famous of them is probably the Barenblat 
one [3]. Thus it is a common practice to use rather rough analytic presentations of 
D{9) to make a qualitative analysis of diffusion process (see, e.g., [3]). 
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In the present paper we propose a new method for description of time evolution of 
a diffusion process and calculation of the diffusion coefficient. The distinct feature of 
our approach is that we find both functions 6 = 6{t, x) and D = D[9) in an explicit 
form, i.e., solve both problems 1) and 2) analytically. To achieve this goal we start 
with experimental data for a particular diffusion system and make the error func- 
tions expansion (ERFEX) of concentration profiles. Analytic description of diffusion 
processes is very convenient for their qualitative analysis. Moreover, our description 
appears to be rather good quantitatively also; its deviation from experimental data 
does not exceed the inaccuracy of measurements. 

We apply this approach to describe in detail the surface diffusion of Li deposited 
on the molybdenum (112) surface which had been previously covered with a submono- 
layer of dysprosium. One more process, the diffusion of Dy adsorbed on Mo(112), is 
used to examine the method generality. Moreover, we believe that it can have a much 
wider application area. 



2 Experimental data and symmetries 

Let us start with experimental data representing surface diffusion of Li on the Mo(112) 
surface precovered with a 0.25 monolayer of dysprosium (below we designate it as 
Dy-Mo(112) surface). The data were obtained in ultra-high vacuum using local mea- 
surements of the work function by a contact potential method (see for details). 

A schematic sketch of the method, termed scanning contact potential microscopy, 
is shown in Fig. [H At the beginning of the experiment, we uniformly covered the 
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Figure 1: Probing the surface distribution of adatoms by scanning contact potential 
microscopy, (a) The initial step-like adatom distribution, (b) The adatom distribu- 
tion after surface diffusion. 

clean surface of a (112) oriented Mo single crystal with dysprosium (its surface density 
amounted to 0.25 of a monoatomic layer) and equilibrated it by annealing at T=1100 
K. The Dy adatoms served thereafter as a controllable admixture which could affect 
the diffusion kinetics of lithium [6]. Then, using a semiplane mask (screen) placed 
between the Li evaporator and the prepared substrate, a half of the crystal surface 
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was covered at room temperature with Li while its another half remained clean of Li 
(Fig. l,a). The surface was scanned with an electron beam formed by an electron gun. 
The movable beam was used to record the distribution of the contact potential (work 
function) over the sample surface by P.A. Anderson's (retarding field) method. In 
separate experiments, the work function change was carefully calibrated with respect 
to the absolute surface concentration of adsorbed atoms (adatoms). The calibration 
data served to convert the work function values to adatom concentrations, n. To 
characterize the relative concentration of adatoms with respect to substrate surface 
atoms, we shall use, as it is conventional in surface science, the term "degree of 
coverage" (or, for short, "coverage"). It is defined as 6' = n/nM-, where um is the 
concentration of the substrate surface atoms {hm = 8.3 ■ lO^^cm^^ for the Mo(112) 
surface). The coverage ^ = 1 is usually termed the geometrical monolayer. 

Under the experimental conditions provided in [5], El [7j, there was neither evap- 
oration of the adsorbate into vacuum nor its diffusion (drain) into volume of the 
Mo substrate. Thus the experimental results relate to a case a "pure" surface diffu- 
sion which could be described by equation ([T]). Notice that the length of the crystal 
sample along the diffusion direction was about 10 mm while the extension of the dif- 
fusion zone in the experiments did not exceed 1.2-1.5 mm. Thus the boundary effects 
connected with the finite size of the sample could be neglected 

The initial Li distribution was step-like shaped. Then upon heating the adsorbate 
profile spreads out due to surface diffusion. Since the edge of the step was oriented 
normally to the atomic channels on the Mo(112) surface, the diffusion proceeded 
quasi-one-dimensionally along the channels, i.e. along the [1, 1, 1] direction [5-7]. 

The experiment consisted of a series of measurements in which we recorded the 
time evolution of the coverage profiles due to diffusion at a constant annealing tem- 
perature. At the beginning of each experiment (t=0), a standard (step-like) initial 
coverage profile of the adsorbate was created and recorded on the crystal kept at room 
temperature, at which the adatom mobility is negligible. This profile is labeled t=0. 
Then the crystal sample was annealed at a fixed temperature. From time to time, the 
annealing was interrupted and the crystal quickly cooled down to room temperature 
to record a new coverage profile arising due to diffusion. After that the annealing was 
continued, and so on. In this way we obtained a series of profiles corresponding to 
different annealing times and a constant annealing temperature. Such experiments 
could be repeated at different annealing temperatures to determine the temperature 
dependence of the diffusion kinetics. 

Measurements were made at times t = 0, 1200, 2100, 3600 and 5400 (seconds) at 
stable temperature T = 600K, experimental error in 6 was AO = 0.003. The recorded 
coverage profiles are presented in Fig. [21 

The initial profile (t = 0) is step-like shaped, but technologically it was impossible 
to form an ideal step. The profiles obtained as a result of surface diffusion have a 
common intersection point at x = 0.88,6' = 0.192 and have rather smooth contours. 
Moreover, a convention is used to set 6' = in the region where its value is below the 
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Figure 2: Coverage profiles of Li adsorbed by Dy-Mo(112) at T=600 K: initial, t = 
(1), and measured at t = 1200 s (2), t = 2100 s (3), t = 3600 s (4), t = 5400 s (5). 
The X- coordinate gives distance in mm. 

measurement accuracy. 

Our task is to give a phenomenological theory of the related diffusion process. 
Abstracting from complicated underlying physical effects which are discussed in pa- 
per [H], we will describe the time evolution of the diffusion process and derive the 
dependence of the diffusion coefficient on concentration. 

To achieve our goal, we will exploit symmetries which form the integral part of 
diffusion processes. A precise analysis of the experimental data makes it possible 
to find a specific symmetry which characterizes coverage profiles. Namely, let us 
fix a profile 6 recorded at time t and consider it as a given function of x. Then 
any other profile 6' measured at time t' can be obtained from 6 using the following 
transformation: 



where xq = 0.88 is coordinate of the point common for all coverage profiles. 

This statement can be verified directly or using computer fits to the experimental 
curves. For example, starting with experimental data for t = 3600 s and applying 
transformation ([2]), one can reproduce profiles for t = 1200 s, 2100 s and 5400 s. 

It should be stressed that the found symmetry is rather exact. As a rule, a 
deviation of curves obtained using transformation ([2]) from experimentally measured 
profiles is within the limits of experimental error, and this deviation decreases with 
growing time. For example, comparing experimental data for coverage profile at 
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t = 5400 s (Fig. [2]) and the corresponding values obtained via change we conclude 
that they are very similar, namely, the differences are below the experimental error. 



3 Time evolution derived from symmetries 

The exact meaning of the symmetries discussed in Section 2 is that functions 6{t,x) 
describing coverage profiles are invariant with respect to the following one-parametric 
group of transformations: 

t^t' = e^^t, x-^x' = xe'' + Xo{l - e"), (3) 

where a is a real parameter. Indeed, solving the first of equations for e" and using 
the second equation we come to relations 

Starting with ([3]) and using tools of the classical group analysis, it is possible 
to describe time evolution of profile 6{t,x). Indeed, the infinitesimal operator of 
transformations ([3]) has the form: 

d ^ d ^ d d d 

^ = + ^TT = + X— - xo — , 
at ox at ox ox 

where rj = |^|a=o and ^ = ^|a=o [S]. The invar iance of coverage profiles with respect 
to transformations ([3]) means that 6{t,x) solves the following equation ^: 

Xe{t,x) = 0. (4) 

It follows from (jlj) that time evolution of coverage profiles 6{t,x) is described by 
the following equation: 

86 X — Xq do do y 86 

di ^ 2t 8^ dt^ 2td^' 

where y = x — xq, 6{t, y) = 6{t, y + xq). As an initial condition we can choose one of 
measured profiles, say that one which corresponds to t = t2 = 1200: 

6ih,y) = 62 = 62iy + xo) (6) 

where 62{x) is a function given numerically in Table 1. Henceforth we omit tilde and 
write 6{t,y) instead of 6{t,y). 

Thus we can describe the time evolution of the coverage profiles without a diffusion 
equation. Solving problem ([5]) with condition ([6]) and using numerical data given in 
the Appendix we can find shapes of these profiles for any time. Some of such profiles 
are presented in Fig. [31 

As we see, the very existence of the symmetries in experimental data makes it 
possible to predict shapes of the coverage profiles which will appear at various times 
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Figure 3: Coverage profiles 6{t,x) for Li on Dy-Mo(112) at T=600 K: experimental 
for t = 2100 s (1) and theoretical for t = 7000 s (2), 10000 s (3), 20000 s (4), 100000 
s (5). The x-coordinate gives distance in mm. 

of heating. This statement is valid for any diffusion system which admits symmetries 
([3]). But if we are interested in analytical description of the diffusion process, we 
have to pose initial conditions analytically. A basic problem is to create a consistent 
model of the diffusion process, i.e. find the dependence of the diffusion coefficient on 
concentration, which is important for physical interpretation. 

In the following sections we solve this problem and give explicit representations 
of coverage profiles 6{t,x) in analytical form. 

4 Symmetries of diffusion models 

Transformations ([2]) have a stable point x = xq = 0.88 and 6 = 0.194 which was the 
only one not being changed. It is convenient to choose just xq as a zero point of our 
coordinate system, i.e. to use variable y — X — Xq — X — 0.88 instead of x. 

Taking into account that in fact we deal with a process which is one-dimensional 
with respect to spatial variables, it is possible to reduce equation ([T]) to the following 
form: 



where D = D{9), y = x — xq. 

Formula ([7]) represents a rather complicated non-linear partial differential equa- 
tion. In addition, we do not know the ^-dependence of the diffusion coefficient D. 
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Fortunately, this equation has a very useful symmetry with respect to scaling of in- 
dependent variables, being invariant with respect to transformations (E]), or 

t^t' = e^^t, y^y' = ye", 6^6' = 6. (8) 

Just this nice property of the diffusion equation makes it possible to use the 

y 

Boltzmann variable ^ = —j^ and search for its similarity solutions 9 = 9{^) where 

both 6 and ^ are invariants of transformations ([H]). 

Notice that equation ([7]) is invariant also with respect to shifts of independent 
variables 

t^t + b, y^y + k (9) 

with arbitrary real parameters h and k. This invariance allows one to choose arbitrary 
initial time and justifies the transition from x io y which we made above. For some 
particular functions D{6), symmetry of the Fick equation is more extended. A com- 
plete classification of symmetry groups of equation has been made by Ovsiannikov 
[To] , a complete group classification of systems of two diffusion equations with source 
terms can be found in [llj and [12j . 

However, symmetries should be compatible also with the initial data of our 
problem, which are of the form: 

e{0,y)=9,{y) , (10) 

where 6i is the initial coverage profile {t = ti = 0) represented numerically in the 
Appendix and graphically in Fig. [21 

We see that the experimentally created profile at t = is not strictly step-shaped, 
and so is not invariant with respect to scaling However, the profile 6i is not far 
from a step and can be considered as a perturbed Heaviside function H{—y) multiplied 
by 0.327 {Omax = 0.327 is the maximum coverage in the initial 6 profile). On the other 
hand, and it is an experimental fact, for sufficiently large times t solutions of our 
problem indeed are invariant with respect to transformations ([2]). And if we apply 
this transformation to infinitely small t', all profiles 9i,92---0i tend to a step-shaped 
one. So we have a direct experimental confirmation of the known mathematical fact 
that similarity solutions of the diffusion equation ([7]) can serve as attractors for other 
solutions. In other words, the yet unknown diffusion coefficient should depend on 9 in 
such a way that the related Cauchy problem ([7]), (fTOl) be asymptotically stable with 
respect to small perturbations of the initial data (see [13] for exact definitions). 

Thus, instead of the actual initial data presented in the Appendix and Fig. [2], we 
can consider an idealized situation when the initial coverage has a step shape, and 
to suppose that for t = the concentration is proportional to the Heaviside function 
H{-y): 

^(0,2/) = e.ax-if(-y) = { (11) 
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where ^max = 0.327 is the maximum concentration in the initial coverage profile. 

The initial- value problem (I7j), (ITTl) is invariant with respect to transformations (IE]), 
and so it is possible to search for invariant solutions 9 = 6{C,) depending on invariant 



Remember that we do not know yet the dependence of the diffusion coefficient D 
on degree of coverage 6. Using experimental data describing dependence of 6' on x at 
fixed heating time t and applying the BM approach [2] it is possible to calculate D 
numerically. Unfortunately, such calculations cannot be done with a sufficiently good 
accuracy, especially in the region of small concentrations. In addition, in this way we 
cannot find the diffusion coefficient in an analytical form. In Section 6 we suggest 
another approach which presupposes direct analytical modelling of profiles 0{^). 

5 Generalized diffusion equation 

Thus we had formulated a possible model of the analyzed diffusion process which is 
based on equation d?]) whose solutions should satisfy the idealized initial conditions 
( TTTl) . However, the possibility of using Fick equation is nothing but a supposition 
which needs additional justifications. In particular, it is necessary to be ensured that 
the diffusion flow is linear in the concentration gradient. 

However, we can use the well justified fact, i.e., the invariance of experimental 
coverage profiles with respect to transformations ([2]), and set the following problem: 
to find the most general evolution equation which is compatible with this symmetry. 

Thus let us suppose that the evolution equation admits symmetries ([2]) and also 
shifts (Q (i.e., does not depend explicitly on space and time variables). Then using 
tools of classical Lie analysis [3], we easily find its general form: 



where P is a dependent variable whose evolution we need to describe, D and G are 
arbitrary functions of P. In particular P can represent a degree of coverage, in this 
case it should be changed to 9. 

We do not present the related routine calculations since our rather strong re- 
strictions (equation should be of evolutionary type and admit the above mentioned 
symmetries) in fact reduce the problem of deducing of ffTSj) to direct use of the di- 
mension analysis. 

Equation ( |T3l) is a direct generalization of the Fick equation (I7j), and the latter 
corresponds to a particular choice G = 0. 





(12) 



^(_oo) = 0.327, 9{+oo) = 0. 




(13) 
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Thus starting with symmetries ([2]), which can be found in the experimental data, 
and symmetries Q, which are natural transformations for the considered diffusion 
system, we deduce the generalized diffusion equation ( fT3l) which includes the Fick 
equation ([7]) as a particular case. Possible physical motivations for generalization of 
equation ([7]) to ( |T3l) are discussed in Section 10. 



6 Modelling functions for coverage profiles 

It is well known that in the case when diffusion coefficient D is independent on 
concentration, the general solution of the problem ([7]), (ITT]) is given by the following 
equation: 

e{t,x) = aerfc(60 

where a = 9^^/2,h = 1/(2V^), i = x/\/t and erfc is the complementary error 
function: 



erf 0(2;) = 1 — erf (2;), erf (z) = 




This fact suggests using just the error function as a constructive element of concen- 
tration profiles for 6'-dependent diffusivities. This idea appears to be very successful 
for description of diffusion processes in general and the processes discussed in Sections 
2 and 9 in particular. 

In this section we present and discuss examples of modelling functions for the cov- 
erage profiles. An algorithm for constructing such functions is given in the following 
section. 

First we consider a rather straightforward representation of profiles 9{S,) which, 
however, is valid only in a reduced interval of the Boltzmann variable ^: 

= %) =ai(l-erf(6iO') (14) 

where oi and hi are parameters. Asking for minimal mean-square deviation of function 
and using MAPLE tools we fix parameters oi and hi to be 

ai = 0.175, 61 = 0.375. 

Function ^(i) perfectly describes the shape of profile 9{^) for ^ lying in the interval 
0.6 < ^ < 00, see Fig. HI In this interval the discrepancy 16* — 6{i)\ does not exceed 
inaccuracy of measurements. However, for ^ < 0.6 this discrepancy increases. 

In order to obtain a modelling function for all non-negative ^, it is sufficient to 
add a small extra term to 6'(i) and define: 

= %)+a2erfc(620, ^>0 (15) 
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Figure 4: Experimental profile ^(^) for t = 5400 (full curve) and the curve 6i given 
by relation (fT^ (broken curve). Units for ^ are 10~^ (mm/s^/^) 

where 02 = 0.02 and 62 = 1-7. 

Function (JTSll gives a very precise presentation for the profile obtained experi- 
mentally at t = 3600 s; the deviation \6 — does not exceed the inaccuracy of 
measurements. Fig. [5] presents the experimental curve and curve defined by equation 
ffT^ . and it is seen that they practically coincide. More exactly, the deviation of 
values of function ffT^ from experimental data is less than 0.003. 
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Figure 5: Coverage profiles 0{C,): experimental data for t = 3600 s (full curve) and 
function given by relation (ITSl) (broken curve). Units for ^ are 10~^ (mm/s^''^) 

One more possible modelling function is given by the following equation: 

^(3) = a[{l - erf (6;0') + a',eTfc{b',0 (16) 

where a[ = 0.145, b[ = 0.395, a'^ = 0.051, b'^ = 0.795. 

The modelling function ffTUj) is a bit more exact but also more complicated. Its 
deviation from the experimental data for t = 5400 s is less than 0.0015, i.e., twice less 
than the experimental error. 
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Thus we wind the modelhng functions fllSI) . 0161] which can be obtained starting 
with a priori supposition that the coverage profiles can be described by second or 
third order polynomials of error functions and demanding minimal root-mean-square 
deviation of these polynomials values from the coverage profiles. This supposition is 
not necessarily valid for other diffusion systems, e.g., it does not give a constructive 
way to build approximate modelling profiles for the system discussed in Section 9. 

In the following section we present a regular way to calculate modelling functions 
for any diffusion system using the error function expansion (ERFEX). 

7 An algorithm for calculation of modelling 
functions 

We have shown above that erfc functions can be successfully used as construction ele- 
ments of the modelling functions of coverage profiles at least for a particular diffusion 
process. Now we shall give a regular way for constructing such functions which has a 
much more extended application area. 

Of course a concrete form of the modelling functions strongly depends on the 
diffusion system, and we cannot propose a universal method how to obtain the most 
simple and exact analytical form of an arbitrary concentration curve. Nevertheless, 
in this section we give an algorithm for calculation of the modelling functions which 
can be applied to any sufficiently smooth concentration profile in the diffusion zone. 

In general case we propose to use ERFEX and search for modelling functions in 
the form 

n 

^(0 = $^A,erfc(fc,(^-e.)) (17) 

i=l 

where C,i {i = 1, 2, ...n) are some fixed values of Boltzmann coordinate, ki and Ai are 
parameters which should be specified. 

Let < or ^1 < ,^2 < • ■ • < ^n, then optimal values of fcj lie inside the interval 

0.25 ,1 

<k,<- -. 18 



6+1 ~ 6 6+1 ~ 6 

In particular it is possible to choose the points 6^ 6 ' ' ' ; Cn in a regular way, i.e., 
with the same distances ^j+i — C,i for all i, and to fix all parameters to be ki = ki = 
■ ■ ■ = kn = p with some p compatible with (1181) . 

Let 6i, 62. ..On be known values of coverage at points 6? ^2? and M be a matrix 
whose elements are Mij = erf c(/cj(^j — ^j)). Then parameters Ai, A2, An are easily 
found by solving the following system of linear algebraic equations: 

M.jAj = ei, z= 1,2, ■■■n, (19) 
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where summing up over the repeated index j is imposed from j = 1 to j = n. 
Equations ( fT9l) are nothing but relation (fTTj) considered at points ^1,^2, ■■■,C,n- 

Let us apply the algorithm to find a modelling function for the diffusion pro- 
cess discussed in Sections 2-6. Consider the results presented in the Appendix for 
t = 5400 s. For simplicity we choose points ^1,^,2, ■■ in a regular way. Namely, we 
chose in the table selected data corresponding to x = 0.63, 0.68, 0.73, 0.78, 0.83, 
0.88, 0.93, 0.98, 1.03 and 1.118 (remember that xq = 0.88), and fix ki = 0.7 for all 
i = l,2,...n (n = 12). At that, system ( fT9l) is easily solved with using symbohc 
calculus program MAPLE. As a result we find the following representation for 6{^): 



e{^) = -0.0086erf c(0.7^ + 3.8865) + 0.0208erf c(0.7^ + 3.2333) 
-O.Olllerf c(0.7^ + 2.5801) + 0.0198erf c(0.7^ + 1.9269) 
+0.0103erf c(0.72; + 1.2737) + 0.0126erf c(0.7x + 0.6205) 
+0.0242erf c(0.7^ - 0.0359) + 0.0328erf c(0.7^ - 0.7544) 
+0.0385erf c(0.7^ - 1.4730) + 0.0211erf c(0.7a; - 2.1915) 
-0.0048erf c(0.7^ - 2.9000) + O.OOlOerf c(0.7^ - 3.6285). 



This a bit cumbersome formula appears to be very precise; in the interval —6 < 
^ < 6 it reproduces experimental data (presented in the Appendix) with an accuracy 
not worse than , moreover, for the majority of experimental points this accuracy is 
not worse than 0.001. Moreover, using all points given in the Appendix for a fixed 
time t, it is possible to find "the most exact" modelling functions 6e whose values 
simply coincide with the experimental curve. However, taking into account the value 
of experimental error, such business does not seem to be reasonable. 

Notice that using the algorithm with a non-regular distribution of points it is 
possible to find a more simple form of the modelling function: 



Ai = 0.103, A2 = 0.234, ^3 = 0.114, = 0.044, 

5i = -1.87, B2 = -0.73, B-i = 1.154, ^4 = 3, R = 0.64. ^ ^ 

Being much more simple than ( |20i) . function (I2T1) is rather exact too; the corre- 
sponding standard quadratic deviation from experimental results is less than 0.003. 

8 Calculation of diffusion coefficient 

Thus we have at our disposal analytic expressions for a close approximation of cov- 
erage profiles. Using them we can calculate the diffusion coefficient D{6) by direct 
integration of equation (fT^ . 



4 




(21) 



i=l 



where 
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First using the simplest modelling function fll4p we find D{6) for concentrations 
9 < 0.16, the related values of ^ satisfy ^ > 0.6. Substituting ffT^ into ffT^ and 
integrating the resultant expression from ^ = z > 1 to infinity, we obtain 

= ( erf (6i^)e"(''i")' + 462erf (6i2)e-(^^")'D) = 0. 

(23) 

Solving ( |23l) for D and using ( 1141) we obtain: 





D = ^|l + ^!^^^i^j, Z = erfinv|,/1--1. (24) 

where erf inv(-) is the inverse error function defined by means of 6' = erf (erf inv(^)). 

Formula (IMl) defines the diffusion coefficient D as an explicit function of concen- 
tration and is valid for all 6 lying in interval < < 0.16. 

In an analogous way, i.e., by direct integration of equation ffT^ . it is possible to 
find the diffusion coefficient D starting with other modelling functions 6{C,) found in 
Sections 6 and 7. The general expression for D is given by the following equation: 



DiO = '\,e% ^ (25) 

which, together with a modelling function for 9{S,)i determines the diffusion coefficient 
as a function of 9 given in a parametric form. 

In contrast to the BM approach, to find the dependence of the diffusion coefficient 
on concentration we simple need to calculate a definite integral of the known function 
and divide it by the redoubled derivative of (known) function with respect to ^. All 
these operations are easily handled using, e.g., MAPLE tools. 

In Fig. we present a plot of the diffusion coefficient fl25p versus the degree of 
coverage 6 given by relation (1201) . These D values are consistent to within 10-25% 
with the data obtained in work [6]. However, the maximum at ^ 0.25 was not 
revealed in [6], where the graphical Matano's evaluation of the coverage profiles was 
applied. This demonstrates that the ERFEX approach provides a more accurate 
processing of experimental diffusion results. 
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Figure 6: Diffusion coefficient D{6) (in units mm^/s ■ 10 ^) of Li on Dy-Mo(112) 
calculated using modelling function (!20ll . T=600 K. 

9 Modelling functions for coverage profiles of Dy 
adsorbed on Mo(112) 

We have seen that modelling functions (l20l) . ( |2T1) give very exact analytic expressions 
for coverage profiles of Li adsorbed on Mo(112). A natural question arises whether it 
is possible to find such functions for modelling the profile shapes for other systems. 

In this section we apply the ERFEX method to another adsorption system, namely 
to Dy adsorbed on Mo(112). Experimental data for this system were obtained and 
discussed in [7] . The related plots of coverage profiles are presented in Fig. [71 

We see that all the profiles recorded in the diffusion process again have a common 
intersection point this time at x = 1.69 which, however, lies out of the initial profile. 
Their shapes are much more specific than ones given in Fig. [21 There are two zones 
with a quick change of coverage and three zones where this change is rather slow. 
These profiles mirror a structural self-organization in the diffusion region, i.e. forma- 
tion of a series of two-dimensional adsorbate phases which differ from each other by 
diffusion parameters and mechanisms ^ . Nevertheless, it appears possible to describe 
these profiles analytically. 

First we represent the coverage profiles using the Boltzmann variable ^ = s/^/t 
and setting the reference frame xq = 0. As a result we conclude that profiles for 
t = 2400 s and t = 4800 s became rather close, the mean quadratic deviation of them 
is less than 0.003. Thus it is possible to describe time evolution of profiles using the 
approach discussed in Section 3. 
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Figure 7: Coverage profile of Dy adsorbed by Mo(112) at T=800 K: initial, t = (1) 
and measured at t = 360 s (2), t = 2400 s (3), t = 4800 s (4) ^. 

Using again the ERFEX, we find the following representation for the coverage 
profile measured at t = 4800 s: 

e{i) = 0.1673erf c(2.2124^ - 6.9912) + 0.06erf c(0.9434^ - 3.3962) 
+0.0094erf c(0.3334^ - 3.0267) + 0.008erf c(4.4444^ - 8.6222) 
+0.006erf c(2.7778^ - 3.4722) + 0.0058erf c(1.21950 + 0.0095erf c(0.817^ (26) 
+2.018) + 0.0102erf c(0.817^ + 5.1716) + 0.0085erf c(1.9608^ + 15.6862) 
+0.062erf c(2.5^ + 22.2) + 0.004erf c(1.9608^ + 19.2157). 

Function fl2U|) is a rather precise approximation of the coverage profile 6{C,) at t=4800 
s; the mean quadratic deviation from the experimental results is less than 0.0015. A 
plot of the related curves is given by Fig. [HI 

We see that ERFEX provides a very good analytical representation for a compli- 
cated coverage profile of Dy adsorbed on Mo(112). 

10 Discussion 

The theory of diffusion is both very old and good developed [H] . Nevertheless, it still 
contains a lot of unsolved problems which attract attention of numerous investigators. 

In the present paper we study three aspects of this theory: using of symmetries 
in experimental data to describe the time evolution of a diffusion process, search- 
ing for a generalized Fick equation which is compatible with these symmetries, and 
construction of modelling functions to describe concentrations of diffusing substances 
and calculate the diffusion coefficient. 

Like the diffusion theory, the basic branch of mathematics which deals with sym- 
metries, i.e., the theory of continuous groups, is rather old too. It was started by 
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Figure 8: Coverage profile of Dy adsorbed on Mo(112) versus the Boltzmann variable, 
measured at t = 4800 s (full curve) and profile described by function ( |26l) (dotted 
curve). T=800 K. 

Sophus Lie around 20 years after appearance of the Fick theory. 

The classical group analysis and its modern generalizations present effective tools 
for investigation and applications of symmetries of mathematical models, including 
diffusion ones. However, to apply these tools directly it is necessary to start with a 
model formulated in terms of (partial or ordinary) differential equations. 

A specificity of diffusion systems is that in general we do not know evolution equa- 
tion a priory. Even if there is a cogent argumentation that the process is described by 
the Fick equation, the dependence of diffusion coefficient on concentration is usually 
unknown. Thus to apply the Lie analysis it is reasonable to search for symmetries in 
experimental data. And it is the first idea which we use in the present paper. 

The second idea is to use erfc functions as a constructive elements of modelling 
functions for concentration curves, i.e., to apply the ERFEX expansion. There are 
two origins of this idea: first, a similarity to solutions of the Fick equation with 
a constant diffusion coefficient, where the erfc function appears very naturally, and 
secondly, the specific shape of this function which seems to be ideal for its use as a 
"brick" for building typical concentration curves. 

Our research is based on a particular diffusion system, i.e., Li adsorbed on Dy- 
Mo(112), which is well studied and described in papers [5]- [8]. This system represents 

^The background of the continuous group theory was formulated by S. Lie in 1875 and pubhshed 
in 1876 [H] 

^Integral, difference and fractional differential equations also are subjects of modern group anal- 
ysis 
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many basic features of surface diffusion, and thus we believe that results of our analysis 
in fact are rather general. A confirmation of this statement was obtained in Section 9, 
where we considered one more surface diffusion system, i.e., Dy adsorbed on Mo(112). 

A precise analysis of the experimental data made it possible to find symmetry 
transformations ([2]) which connect experimental curves giving the coverage profiles. 
Then we apply these symmetries to describe time evolution of the diffusion system 
and to formulate a possible generalization of the Fick equation, given by relation ( fT3l) . 
which keeps invariance with respect to transformations ([3]). A physical motivation 
to search for such generahzations is caused by the fact that the values, measured 
immediately in the experiment to judge of the concentration of diffusing species, are 
not necessarily in direct proportion to the concentration. Moreover, in some cases we 
may actually be interested not just in the concentration of particles, but rather in 
various physical properties connected with it (e.g. electrical conductivity , mechanical 
strength, optical properties, work function etc.). Suppose that concentration 6 and 
the relevant system property P are related by the dependence 

e = F{P). (27) 

Substituting (1271) into the Fick equation ([7]) we come to equation ( |T3l) . where 

D{P)=D{F{P)), G = D^, ^' = W ^^^^ 

On the other hand, if a diffusion process is accompanied by dissipation (e.g. due to 
evaporation or chemical reaction of the diffusing species), the related mathematical 
model also needs a generalization of the Fick equation by inclusion of the terms 
do 

depending on — . In order to create a mathematical model of such process which 
ox 

keeps symmetries (E]), one should use just the motion equation ( fT3l) with P = 6. 

We see that the generalized equation (fT3l) appears as a result of rather straightfor- 
ward considerations and presents an alternative to the generally recognized equation 
([7]). Surely, for the case when ^ is a linear function of P and there is no dissipation, 
equations ( fT3l) . ( !28|) are reduced to ([7]). 

Use of modelling function for coverage profiles is very convenient for analysis of 
diffusion processes. These functions are sufficiently exact and can be used for direct 
calculation of diffusion coefficients starting with the Fick equation. 

In the present paper we propose a few modelling functions for coverage profiles 
of Li adsorbed on Dy-Mo(112). They described experimental data with a high preci- 
sion, the deviation from experimental data is less than experimental error. Thus the 
modelling functions present useful tools for both qualitative and quantitative studies 
of the diffusion systems. 

Let us note that calculating the modelling function for small concentrations with 
higher accuracy (i.e., using all experimental given points) it is possible to find a spike 
for the diffusion coefficient for small concentrations. The related plot is given by Fig. 
9. 
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Figure 9: Diffusion coefficient D{9) (in units inmVs • lO"*^) of Li on Dy-Mo(112) 
calculated precisely for small concentrations. T=600 K. 

Notice that the spike in Fig. 9 is indicated for a very small concentrations com- 
patible with the experimental error and so it cannot be treated as well experimentally 
justified. On the other hand this spike can be observed for all series of experimental 
data presented by Fig. 2 and also in some other diffusion systems, e.g., in Dy ab- 
sorbed on Mo(112). Indeed the coverage profile given by Fig. 8 is rather steady at 
^ = 5 — 7 when 6 ~ 0.02 and so the related diffusion coefficient will have a maximum 
thanks to small value of the derivative dO/d^. 

Concerning interpretation the spike indicated by Figs. 9 and min.-max. presented 
by Fig. 6 we can mention that the coverage dependence of diffusivity is due to lateral 
interactions of diffusing atoms and partially also to the specific features of the sub- 
strate atomic structure (both intrinsic and caused by defects). The combined action 
of lateral interactions and surface potential corrugation determines the sequence of 
phase transitions that occur in the diffusion zone. The phases can differ from one an- 
other not only by the diffusion parameters, but also by the diffusion mechanisms (see 
e.g. a review [8] and references therein). The maxima of diffusivity at low coverage 
were observed for a number of systems. We attribute this effect to a high mobility 
of single adatoms. The fast decrease of D with growing coverage may be caused 
by formation of clusters, whose diffusion mechanisms can be very diversified, but 
their mobility is generally lower than that of single adatoms. The diffusivity is also 
rather low in the regions of first-order phase transitions in adlayers. As the coverage 
grows further and the adlayer becomes increasingly dense approaching a close-packed 
monolayer, the diffusion takes on a pronounced collective character. In particular. 
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in the region of commensurate-incommensurate phase transition the diffusion seems 
to proceed by a relay-race walks of misfit dislocations (topological solitons), which 
provides a high diffusion rate (a local D maximum). This effect was revealed for many 
adlayers. For more details refer to papers [8], [T6] . 

Summarizing, we propose a constructive and convenient algorithm (ERFEX) for 
generating of modelling functions which is valid for arbitrary sufficiently smooth 
curves not necessarily related to a diffusion process. In particular, using this al- 
gorithm and starting with experimental data, it is possible to determine the diffusion 
coefficient with a higher accuracy than the BM approach and spline approximation. 
The algorithm can be treated as a specific generalization of the wavelet approach 
which can be applied to study of diffusions. 
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Appendix 

Here we present the table including experimental data obtained for coverage pro- 
files of Li adsorbed by Dy-Mo(112). They are used in the main text to estimate 
exactness of the modelling functions. 
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Coverage profiles of Li adsorbed on Dy-Mo(112). 



X (mm) 


6'i,t = 


92,t= 1200 s 


6^3, t = 2100 s 


6*4, t = 3600 s 


6*5, t = 5400 s 





0.327 


0.327 


0.327 


0.327 


0.327 


0.20 


0.323 


0.323 


0.324 


0.322 


0.322 


0.40 


0.322 


0.323 


0.319 


0.318 


0.316 


0.48 










0.310 


0.56 










0.300 


0.60 


0.317 


0.316 


0.309 


0.300 


0.291 


0.64 




0.310 


0.305 


0.293 


0.285 


0.68 


0.317 


0.317 


0.295 


0.283 


0.274 


0.72 


0.316 


0.297 


0.288 


0.272 


0.262 


0.76 


0.318 


0.283 


0.270 


0.258 


0.248 


0.80 


0.312 


0.266 


0.253 


0.239 


0.234 


0.84 


0.299 


0.232 


0.221 


0.217 


0.215 


0.88 


0.201 


0.193 


0.194 


0.194 


0.195 


0.90 


0.067 


0.165 




0.177 


0.181 


0.92 


0.036 


0.137 


0.155 


0.159 


0.169 


0.94 


0.025 


0.103 


0.131 


0.144 


0.154 


0.96 


0.013 


0.066 


0.104 


0.122 


0.140 


0.98 


0.010 


0.040 


0.074 


0.105 


0.124 


1.00 


0.008 


0.018 


0.051 


0.081 


0.106 


1.02 




0.012 


0.030 


0.067 


0.091 


1.04 


0.003 


0.008 


0.018 


0.049 


0.077 


1.08 


0.001 


0.003 


0.005 


0.021 


0.047 


1.12 





0.002 


0.003 


0.009 


0.027 


1.16 












0.011 


1.2 











0.001 


0.004 


1.28 










0.001 
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